Use marmap::getNOAA.bathy() to load in New Zealand bathymetry data directly from NOAA.
nz <- getNOAA.bathy(
lon1 = 162,
lon2 = 180,
lat1 = -33,
lat2 = -50,
resolution = 1
)
## Querying NOAA database ...
## This may take seconds to minutes, depending on grid size
## Building bathy matrix ...
Use sf::read_sf() to read the shapefile containing the New Zealand coastline.
nz_coast <- sf::read_sf(
here(
"data",
"nz-coastlines-and-islands-polygons",
"nz-coastlines-and-islands-polygons-topo-150k.shp"
)
)
Base R plot
blues <- c("lightsteelblue4", "lightsteelblue3",
"lightsteelblue2", "lightsteelblue1")
greys <- c(grey(0.6), grey(0.93), grey(0.99))
plot(nz, image = TRUE, land = TRUE, lwd = 0.03,
bpal = list(c(0, max(nz), greys),
c(min(nz), 0, blues)))
# Add coastline
plot(nz, n = 1, lwd = 0.4, add = TRUE)
Using the marmap::autoplot.bathy() function to plot with ggplot.
autoplot.bathy(nz, geom=c("tile","contour")) +
scale_fill_gradient2(low="dodgerblue4", mid="gainsboro", high="darkgreen") +
labs(y = "Latitude", x = "Longitude", fill = "Elevation") +
coord_cartesian(expand = 0)+
ggtitle("A marmap map with ggplot2")
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Convert raster object to dataframe for ggplot.
nz %>%
as.raster() %>%
rasterToPoints() %>%
as.data.frame() ->
nz_df
Plot New Zealand Coastline with bathymetric and elevation data.
ggplot() +
geom_raster(data = nz_df
# ## Uncomment below to only get bathymetry
# %>%
# filter(layer <= 0)
,
aes(x = x, y = y, fill = layer)) +
# ## I think it looks better without the contours on land, if any at all
# ## Uncomment first and last two lines below to get all contours
# ## Uncomment all below lines to get only seafloor contours
# geom_contour(data = nz_df
# %>%
# filter(layer <= 0)
# ,
# aes(x = x, y = y, z = layer), color = "grey20", size = .1) +
scale_fill_viridis_c(option = "A") +
geom_sf(data = nz_coast, fill = NA, color = "black", size = .15) +
coord_sf(xlim = c(162, 180), ylim = c(-33, -50), expand = c(0, 0)) +
labs(fill = "Elevation\n&\nBathymetry",
x = "Longitude",
y = "Latitude",
title = "New Zealand",
subtitle = "Digital elevation (1 arc minute resolution)") +
theme_classic()
Convert from bathy object to RasterLayer Object
nz %>%
as.raster() ->
nz_raster
Plot Raster onto interactive leaflet map
leaflet() %>%
addProviderTiles(providers$Esri.OceanBasemap, group = "Ocean Base Map") %>%
addRasterImage(nz_raster, group = "Bathymetry") %>%
addPolygons(data = nz_coast, weight = 1, opacity = 1, fill = FALSE) %>%
addLayersControl(
baseGroups = c("Ocean Base Map"),
overlayGroups = c("Bathymetry"),
options = layersControlOptions(collapsed = FALSE)
)